# Replication Archive for 
# Alexander Coppock, Emily Ekins and David Kirby (2018), 
# "The Long-lasting Effects of Newspaper Op-Eds on Public Opinion",
# Quarterly Journal of Political Science: Vol. 13: No. 1, pp 59-87.

rm(list = ls())

# Uncomment to install
#install.packages(c("dplyr", "stargazer", "MASS", "estimatr"))

library(dplyr)
library(stargazer)
library(estimatr)
# MASS not loaded but called for MASS::polr()

# Set working directory to replication archive

load("mturk_opeds_cleaned.rdata")
load("elite_opeds_cleaned.rdata")
load("dv_df.rdata")
source("opeds_source.R")

fit_1 <- lm(dv_amtrak_main_w1 ~ Z, data = filter(mturk_opeds, responded_w1))
fit_2 <- lm(dv_climate_main_w1 ~ Z, data = filter(mturk_opeds, responded_w1))
fit_3 <- lm(dv_flat_main_w1 ~ Z, data = filter(mturk_opeds, responded_w1))
fit_4 <- lm(dv_vets_main_w1 ~ Z, data = filter(mturk_opeds, responded_w1))
fit_5 <- lm(dv_wall_main_w1 ~ Z, data = filter(mturk_opeds, responded_w1))

fit_6 <- lm(dv_amtrak_scale_w1 ~ Z, data = filter(mturk_opeds, responded_w1))
fit_7 <- lm(dv_climate_scale_w1 ~ Z, data = filter(mturk_opeds, responded_w1))
fit_8 <- lm(dv_flat_scale_w1 ~ Z, data = filter(mturk_opeds, responded_w1))
fit_9 <- lm(dv_vets_scale_w1 ~ Z, data = filter(mturk_opeds, responded_w1))
fit_10 <- lm(dv_wall_scale_w1 ~ Z, data = filter(mturk_opeds, responded_w1))

stargazer(fit_1,fit_2, fit_3, fit_4, fit_5,
          se = starprep(fit_1,fit_2, fit_3, fit_4, fit_5),
          p = starprep(fit_1,fit_2, fit_3, fit_4, fit_5, stat = "p.value"),
          style = "apsr", column.sep.width = "0pt",  digits = 3, 
          omit.stat = c("adj.rsq", "f", "ser"), 
          dep.var.labels = c("Amtrak", "Climate", "Flat Tax", "Veterans", "Wall Street"), 
          covariate.labels = c(
            "Op-ed: Amtrak",
            "Op-ed: Climate",
            "Op-ed: Flat Tax",
            "Op-ed: Veterans",
            "Op-ed: Wall Street",
            "Constant (Constant)"
          ), 
          model.numbers = FALSE, 
          notes = c("Models estimated via OLS. Robust standard errors in parentheses."), 
          notes.append = TRUE, font.size = "small", float = TRUE, align = TRUE, label = "tab:mturk_main", 
          title = "MTurk Experiment: Treatment Effects on Main Dependent Variables")

stargazer(fit_6, fit_7, fit_8, fit_9, fit_10,
          se = starprep(fit_6, fit_7, fit_8, fit_9, fit_10),
          p = starprep(fit_6, fit_7, fit_8, fit_9, fit_10, stat = "p.value"),
          style = "apsr", column.sep.width = "0pt",  digits = 3, 
          omit.stat = c("adj.rsq", "f", "ser"), 
          dep.var.labels = c("Amtrak", "Climate", "Flat Tax", "Veterans", "Wall Street"), 
          covariate.labels = c(
            "Op-ed: Amtrak",
            "Op-ed: Climate",
            "Op-ed: Flat Tax",
            "Op-ed: Veterans",
            "Op-ed: Wall Street",
            "Constant (Constant)"
          ), 
          model.numbers = FALSE,
          notes = c(
            "Models estimated via OLS. Robust standard errors in parentheses.",
            "Dependent variables are constructed via factor analysis and have unit variance."
          ), 
          notes.append = TRUE, font.size = "small", float = TRUE, align = TRUE, label = "tab:mturk_scale", 
          title = "MTurk Experiment: Treatment Effects on Composite Scale Dependent Variables")

# Elite Experiment --------------------------------------------------------

fit_1 <- lm(dv_amtrak_main_w1 ~ Z, data = filter(elite_opeds, responded_w1))
fit_2 <- lm(dv_flat_main_w1 ~ Z, data = filter(elite_opeds, responded_w1))
fit_3 <- lm(dv_vets_main_w1 ~ Z, data = filter(elite_opeds, responded_w1))
fit_4 <- lm(dv_wall_main_w1 ~ Z, data = filter(elite_opeds, responded_w1))

fit_5 <- lm(dv_amtrak_scale_w1 ~ Z, data = filter(elite_opeds, responded_w1))
fit_6 <- lm(dv_flat_scale_w1 ~ Z, data = filter(elite_opeds, responded_w1))
fit_7 <- lm(dv_vets_scale_w1 ~ Z, data = filter(elite_opeds, responded_w1))
fit_8 <- lm(dv_wall_scale_w1 ~ Z, data = filter(elite_opeds, responded_w1))

stargazer(fit_1,fit_2, fit_3, fit_4, 
          se = starprep(fit_1, fit_2, fit_3, fit_4),
          p = starprep(fit_1, fit_2, fit_3, fit_4, stat = "p.value"),
          style = "apsr", column.sep.width = "0pt",  digits = 3,
          omit.stat = c("adj.rsq", "f", "ser"),
          dep.var.labels = c("Amtrak", "Flat Tax", "Veterans", "Wall Street"),
          covariate.labels = c("Op-ed: Amtrak", "Op-ed: Flat Tax", "Op-ed: Veterans","Op-ed: Wall Street", "Constant (Constant)" ),
          model.numbers = FALSE,
          notes = c("Models estimated via OLS. Robust standard errors in parentheses."),
          notes.append = TRUE, font.size = "small", float = TRUE, align = TRUE, label = "tab:elite_main",
          title = "Elite Experiment: Treatment Effects on Main Dependent Variables")

stargazer(fit_5, fit_6, fit_7, fit_8,
          se = starprep(fit_5, fit_6, fit_7, fit_8),
          p = starprep(fit_5, fit_6, fit_7, fit_8),
          style = "apsr", column.sep.width = "0pt",  digits = 3,
          omit.stat = c("adj.rsq", "f", "ser"),
          dep.var.labels = c("Amtrak", "Flat Tax", "Veterans", "Wall Street"),
          covariate.labels = c("Op-ed: Amtrak",  "Op-ed: Flat Tax", "Op-ed: Veterans","Op-ed: Wall Street", "Constant (Constant)" ),
          model.numbers = FALSE,
          notes = c(
            "Models estimated via OLS. Robust standard errors in parentheses.",
            "Dependent variables are constructed via factor analysis and have unit variance."
          ),
          notes.append = TRUE, font.size = "small", float = TRUE, align = TRUE, label = "tab:elite_scale",
          title = "Elite Experiment: Treatment Effects on Composite Scale Dependent Variables")

# Ordered models for R1 ---------------------------------------------------

fit_1 <- MASS::polr(factor(dv_amtrak_main_w1) ~ Z, data = filter(mturk_opeds, responded_w1))
fit_2 <- MASS::polr(factor(dv_climate_main_w1) ~ Z, data = filter(mturk_opeds, responded_w1))
fit_3 <- MASS::polr(factor(dv_flat_main_w1) ~ Z, data = filter(mturk_opeds, responded_w1))
fit_4 <- MASS::polr(factor(dv_vets_main_w1) ~ Z, data = filter(mturk_opeds, responded_w1))
fit_5 <- MASS::polr(factor(dv_wall_main_w1) ~ Z, data = filter(mturk_opeds, responded_w1))

fit_6 <- MASS::polr(factor(dv_amtrak_main_w1) ~ Z, data = filter(elite_opeds, responded_w1))
fit_7 <- MASS::polr(factor(dv_flat_main_w1) ~ Z, data = filter(elite_opeds, responded_w1))
fit_8 <- MASS::polr(factor(dv_vets_main_w1) ~ Z, data = filter(elite_opeds, responded_w1))
fit_9 <- MASS::polr(factor(dv_wall_main_w1) ~ Z, data = filter(elite_opeds, responded_w1))
